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ABSTRACT 

Master's Thesis which discusses nature of allocation 
problems. Danskin Algorithm for solution of a convex func 
tion to be minimized over a closed convex, set is developed 
An example of application involving solution of a 3600 
variable allocation problem using a computer is provided. 
Paper includes analysis of solution and discussion of pro- 
blems encountered in computer application.. 
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I. INTRODUCTION 



A large selection of problems faced within, the mili- 
tary and commercial environments consist of finding the 
"best” way of using available resources. In most cases 
the "best" method of using these resources is that which 
provides the most desirous return to the using agency or 
organization. Typically, this return can be represented 
as a payoff associated with the use- of resources',, or in 
mathematical terms it is a function - the value of which 
is dependent on and determined by the amount or value of 
resources expended. Naturally, the availability and use 
of these resources are bounded and, therefore, act as a 
constraint on the using organization. This constraint 
could be that the organization is required to use some set 
minimum amount of a particular resource, or in the more 
usual sense only a given maximum ai.iount of the resource 
is available for use. In the former case the resource is 
termed lov;er bounded, and in the latter case the resource 
is considered upper bounded. In many such problems the 
resource is both upper and lower bounded and usually some 
allocation of the resource between the bounds is acceptable. 

If the situation exists where the bounds on a resource 
are known but where varying amounts of the resource can be 
expended in a number of different ways, then the manner in 
which the resource is allocated becomes of critical importance 
to the allocating organization in terms of the payoff it 
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receives from choosing a particular allocation scheme. 

By simple use of a crystal ball, random selection device, 
or some other means; and by exercising care not to violate 
the bounds, the organization can adequately allocate. How- 
ever, it cannot be sure it has found the best allocation ' 
scheme . 

If the payoff realized can be expressed as a mathe- 
matical function of the resources allocated, then the pro- 
blem becomes susceptible to solution by matkematical 
programming methods. In many cases then the "best" or 
optimal allocation of resources can be found and proven to 
be optimal. These type problems can be described as "al- 
location" type problems. Typically, the value representa- 
tion of a resource variable is a non-negative number less 
than one and indicating the percent of the total resource 
available which is to be expended in the particular manner 
or activity associated with the given variable. 

The formulation of the payoff function is usually the 
deciding factor in determining what type of mathematical 
programming method could be used to find the optimal al- 
location or optimal solution to the problem. For instance, 
if the payoff function is linear in nature, then it may be 
possible to solve the problem by using existing linear pro- 
gramming methods. Hov^:ever, many payoff functions cannot be 
expressed as linear relationships even when attempts are 
made at reasonable approximations. In this case one must 
turn to other methods. In some cases it may be true that 
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no known mathematical method for optimizing exists. In this 
case it may be proper to reformulate the problem or possibly 
return to the crystal ball. 

This paper will address a particular category of al- 
location problems, develop a method for solving problems 
which fit into this category, and provide a detailed exam- 
ple of a p ractical problem solved by the method developed. 
Credit for development of this method belongs to Dr. John 
M. Danskin, Professor - U. S. Naval Postgraduate School. 

This method, employing an algorithm termed the Direction 
Finding Algorithm, was presented by Dr. Danskin in various 
classroom lectures presented during his tenure at that 
school. 

A. PRESENTATION OF THE PROBLEM 

Consider an allocation type problem in which the pay- 
off function is convex (or concave) and it is decided that 
the optimal solution is that solution which minimizes (or 
maximizes) the payoff function. Now, consider possible al- 
location of resources to be represented as X = (xi ,X 2 , . . .x^) 
where = percent of resources expended in activity #1, 

X 2 = percent of resources expended in activity #2, etc.; for 

X — * 1,..., n . 

As indicated earlier, certain constraints must be heeded 
in finding the optimal solution. These constraints can be 
formulated as follows: 

(1) It is necessary and desirable to use up all of the 
resources available. 
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(2) It may or may not be necessary to use some per- 
cent of the resource in any given activity, 

(3) The percent of resource expended in activity i 
cannot exceed its upper or lower bounds, and 

(4) Negative use of the resource is meaningless and 
not possible. 

This problem as described can then be mathematically 
expressed as: 

maximize f(x) where f(x) is concave 
or 

minimize f(x) where f(x) is convex 

subject to constraints: 

E X. = 1 
1 

1 

a. < X. < b . 

1 - 1 - 1 

where a^ and b^ represent the upper and lower bounds re- 
spectively of the allocation to the i^^ activity. It then 
follows that the following relationships are true: 

0 < a. < b- 

1 1 

E a. < 1 < Z b. . 

1 . 1 

1 1 

Essential to the development and use of the Direction 
Finding Algorithm is the fact that the payoff or objective 
function must be differentiable. Therefore, this require- 
ment that the gradient exists and can be calculated is con- 
sidered as an added constraint. 
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This category o£ allocation problems just described is 
the category for which the about- to-be-developed Direction 
Finding Algorithm can be used to find the optimal solution. 
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II. DERIVATION AND DEVELOPMENT OF THE ALGORITHM 



In this section the Direction Finding Algorithm will 
be developed and a step-by-step procedure described for 
its application. 

Recall the original problem is to maximize (minimize) 
a concave (convex) differentiable function of the vector 
X = (xi ,X 2 , . . . ,x^) with constraints as follows: 

Ex.=l 0<a. <b. 

1 - 1 1 



a. 

1 



< X . < b . 

1 - 1 



E 

i 



a. < 1 < E b . 
1 - 1 
1 



The basis of the algorithm is to find the direction 
of fastest increase. This means it is necessary to maxi- 
mize the Directional Derivative which is defined in Ref. [I] 
as 

D f(x) = E Y- f (x) = E ^-y- 
Y -1 x.^-^ . I'l 

' 111 



where represents the gradient and y represents the direc- 
tion of the gradient. 

To develop the algorithm, the following constraints are 
imposed on y, the direction vector: 



Y^ > 0 
Yi < 0 



E Y- = 0 
. ' 1 
1 

E Y? = 1 
i ^ 

if X. = a. 
1 1 

if X . = b . 
1 1 
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The entire problem then becomes 



maximize 


zn • Y 








subject to 


X Y? 


= 


1 






X Yi 


= 


0 






’fi 


> 


0 


if X. = a. 
1 1 






< 


0 


if X. = b . 
1 1 



with the assumption that the maximum is non-negative. 

The first step to solving this problem will be to show 
necessary conditions for solution by application of the 
well-known Kuhn-Tucker (K-T) Conditions which can be found 
in Ref. [2] and many other places in recent literature on 
non-linear programming. The constraints can be rewritten 
and LaGrange Multipliers assigned as follows: 



CONSTRAINT 




MULTIPLIER 


Z Y? - 1 


< 


0 


X! 

1 


1 - Z Y* 


< 


0 


X! ' 
1 


E Yi 


< 


0 




-2 Yi 


< 


0 


u:- 


-^i 


< 


0 if X. - a. 
1 1 


V . 

1 


Yi 


< 


0 if X- = b . 

1 1 


TT . 
1 



Application of the Kuhn-Tucker Theorem then results in 
V (fi*Y) = X’V(Z y|-1) + X”V(1-Zy|) 

+ p'V(Z Y^) + y”V(-Z Y^) 

+ E'(-v.)V(yP + 2”(Tr^)V(Y^) 
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where 

V as usual represents the gradient 

Z' implies summation over those i's for which 

X. = a. 

1 1 

Z’’ implies summation over those i's for which 



If at this point only the i^h component of the gradient 
is considered, the above equation leads to 



J2. 

1 



Y?(2y[ + 2r[') + 




V . + IT . 

r 1 



where again enters consideration only if = a^ and 
is considered only if x^ = b^. The y^ represents the i^^ 
component of the y vector in the optimal solution as pro- 
vided for in the Kuhn-Tucker Theorem. Letting X = 2X' - 2X ' ' 
and y = y ' - y'', the following results can readily be de- 
duced from the previous equation and the original con- 
straints : 




or if X. = a. and y- > 0 
r 1 1 

or if X. = b. and y- < 0 
1 1 ' 1 

if X. = a. and y- = 0 
if x^ = b^ and Y^ = 0 



Since v. >0 and u. >0 from the Kuhn.-Tucker Theorem, 
1 - 1 - 

it is equivalently true that 
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0 



( 1 ) 



< X Y? + y 

fl. > X Y- U 

1 - 1 

. = X Y • y 

1 ' 1 



if X. = b. and y* = 

1 1 1 

if and Y? - 0 (2) 



otherwise conditions satis- 
fying the original con- (3) 
straints . 



Therefore, if the optimal solution exists, then X and 
y exist and are non-negative and satisfy the above condi- 
tions . 

The next step is to consider the sufficiency of the 
conditions. It is asserted that if X and y exist and X is 
non-negative and if y° satisfies results in (1) , (2) ,. and 
(3) above, then y° maximizes Zn*Y- 

To prove this let y be any arbitrary direction vector 
satisfying the conditions derived via the Kuhn-Tucker 
Theorem. Observe that n*Y = " u)Yj^ since 

= 0 from the original problem constraints. This is true 
for all i representing elements of the n and y vectors. 
Therefore, the summation can be broken into three groups 
of terms as indicated in (1), (2), and (3) above; i ..e . , 



n-Y = E^(fi.-y)Y- + Z^^(fi.-y)Y. + 2^^^(n.-y)Y. 



where and are summations over terms in (1), (2), 

and (3) respectively. < 0 since (Q^-y) < 0 from the K-T 

necessary condition (1) and Yj^ > 0 from the original con- 
straint. By similar reasoning it can be seen that E^^ < 0. 
Therefore, it is true that < 2^^^. 



However, E^^^(fl^-y)Y^ = XE^^^Y^Y^ from K-T condition 

(3) and XE^^^y-Y- ~ ^ in Y-Y- since y- 0 for all i not 

^ 'I'l all'i'i 1 



1 
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contained in T. 



Ill j-v 1“ 0 ^ \ f 02^/'I■ 2 

■ ^ all 5 ^tall ’(all - X 

1 1 ^ i " 

from the Schwartz inequality and from the original con- 
straints. Therefore, fi*y < but = E(f2^-VL)y^ = AZy?^ 

= A. And since A > 0, it is true that fi*y < f2-y° with the 
resultant fact that y® maximizes fi*y as was to be shown. 
Furthermore, the maximum value is A or in terms of the suf- 
ficiency conditions; if A can be found then the optimal 
solution is obtained. 

At this point it is convenient to rewrite equation 
(3) and the conditions necessary for the equality: 



n. = Ay? + y 
1 ' 1 



if a. < X- < b. 
Ill 



or if X. = a. and y? > 0 
1 1 ' 1 



or if X. = b. and y? < 0 
1 1 ' r 



Letting Z' indicate a summation over the terms ful- 
filling these conditions and summing, results in 



or 



Z'fi. = AZ'y? + Z'y = Z'y 
1 ' 1 



y , Z ’ 



(4) 



where N' is the total number of terms fulfilling the condi- 
tions. If the equation is squared first and then summed, 
the expression 



A = Z'(J^. -y)2 

results. It also follows from substitution that 



(5) 



y? = (fi. - y)/A, 



( 6 ) 
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With a method of calculating values for y, X, and 
now in hand, development of the algorithm can proceed. 

For notational ease, let DS = the ’’Distinguished Set” 
and be defined as that set of indices i such that a.<x.<b., 
or X. = a. and y? > 0, or x. = b. and y? < 0. In. other 
words, the distinguished set represents those elements of 
the allocation vector, X, where the element does not lie on 
the boundary; or if the element is on the boundary, then 
the corresponding direction vector • element implies a move 
away from the boundary. With this notation equa.tion (4) 
can be rewritten as 



V 




E 

ieDS 



Q. 



1 



(7) 



where total number of elements in DS .. 

The elements contained in DS can be separated into 
three categories as follows: 

Category #1. The X element is not on its boundary and 
there is no restriction on the y element. 

Category #2. The X element is on its lower boundary 
and the y element is positive. From equation (6) then 

. > y . 

1 

Category #3. The X element is on its upper boundary 
and the y element is negative. From equation (6) then 

fi. < y . 

1 

It is evident now that the algorithm for finding y° , 
or the direction of fastest increase, must search out and 
find each of the elements included in the three categories 
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above. Recall that y is actually the average of all for 
which / 0* This means it is necessary to determine y 
concurrently with determining the elements in DS . The fol- 
lowing procedure provides the method for doing just this. 

A. THE DIRECTION FINDING ALGORITHM 

1 y. 

Step #1. Let y = rr — . ^2. and DS = {<})} to start. 

^DS 1 

Step #2. Place into DS each index i for which a.<x.<b.. 

Note that this implies no restriction on Note also that 

throughout the algorithm, once i has been entered into DS , 

it is no longer considered as a candidate for entry. 

Step #3. If DS = {(})}, then temporarily let y = min{J^.}. 

i 

Do not place this i into DS (the minimum gradient element 
serves simply as a starting point). If DS ^ {<})}, then cal- 
culate y from equation (7) . 

Step #4. Search indices for cases where x^ = b^. Find 

the smallest associated with this search. If this < y, 
1 1 

then enter this i into DS and recalculate y. If this ft > U» 
then proceed to Step #6. 

Step #5. Repeat Step #4 until no additional i's can 
be entered into DS . 

Step #6. Search indices for cases where x^ = a^ . Find 
the largest ft^ associated with this search. If this ft^ > y , 
then enter this i into DS and recalculate y. If this ft^ < y , 
then proceed to Step #8. 

Step #7. Repeat Step #6 until no additional i's can 
be entered into DS. 
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Step #8. Repeat Step #4 through Step #7 until no ad- 
ditional i's can be entered into DS . 

Step #9. Calculate X = [ Z If A = 0 

icDS ^ 

then STOP . 

Step #10. Calculate from equation (6). 

This concludes the algorithm. At this point two 
special case considerations should be discussed; 

CASE I. The possibility exists that only one element 
will be entered into the DS . In this event A = 0 and the 
optimal solution lies on the boundary. y° at this point 
indicates that any better solution would violate the orig- 
inal constraints. 

CASE II. If more than one element has been entered 
into DS and A = 0, then fi. = p for all i. In this event 
all are obviously equal and the optimal solution has 
been reached. In fact, the value of A is a measure of how 
close the procedure has come to the optimal solution. It 
may be true in a practical application of the algorithm that 
A will diminish more slowly with each successive computation 
as the optimal solution in neared. It may be necessary in 
this case for the user to make a determination of the pre- 
cision desired and to terminate calculations when a satis- 
factory degree of precision has been reached. 

In the next section an example application of the algo- 
rithm in a practically oriented problem is provided. 
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III. EXAMPLE APPLICATION OF DIRECTION FINDING ALGORITHM 



Consider a military oriented problem in. which an ord- 
nance delivering organization has responsibility for and 
the capability of firing its ordnance into a particular 
sector of territory. Assume that this organization has 
some means of assigning probabilities to the potential 
existence of a target at any of a number of specific loca- 
tions within its sector of responsibility. Assume also 
that the organization knows what damage it can expect from 
a specific type round delivered at a specific point. Know- 
ing this information, the organization would then like to 
know how to optimally distribute the inventory it intends 
to deliver so that the probability of target survival is 
minimized. ^ 

As a first step to solution of this problem, the sector 
of responsibility is assumed to be square in geometric shape 
and an imaginary grid is imposed upon the sector. The in- 
dividual squares of the grid are numbered sequentially be- 
ginning in the top left corner and working first across and 
then down. Location of a target (or the possible location 
of a target) can then be referenced by a specific numbered 
square of the grid. Expected damage in any grid square can 
be calculated from the number of rounds which impact in that 



1 

This problem is formulated and presented for solution 
in a thesis paper prepared by Captain William A. Hesser, 
USMC, at the U. S. Naval Postgraduate School in February/ 
March 1971. 
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square and/or in nearby squares. Allocation of ordnance 
is represented by the number of rounds delivered to a 
specific grid square. 

The survival probability of a target, or the objective 
function which is to be minimized, can now be expressed as 
a mathematical equation 



F(Yj = Z p exp{-Z B..y.} 
i ] ^ 

where the variables have interpretations as fallows: 

p^ = probability that target is located in i^h square, 

= probability of destruction of a target appearing 
in the i^ square by a round detonating in the 



:th 



square (this is referred to as the damage 



Y: 



function) , 

= percent of total inventory to be expended allo- 

j 

cated to the jth square, 

F(Y) = survival probability of the target and is a 
function of the allocation of Y. 

Note at this point that the survival function, F (Y) , 
is a convex function and the variable Y is actually a 
closed convex set (there are only a finite number of rounds 
which can be fired) . Also note that in the trivial case 
when Y = (0,0,...,0), i.e., no firing takes place, the 
survival function reduces to F (Y) “ ~ ^ which is ap- 

propriate . 

The problem now is to determine what allocation of Y 
will minimize F(Y) . 
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It can be seen at this point that in order to obtain 



realistic and usable results, it is necessary to use a grid 
with a large number of squares. Correspondingly,, the allo- 
cation vector is composed of a large number of elements. 

The problem can now be represented as follows: 

minimize FfY) = Z p . exp { - Eg . .y . } 

j ij 3 

subject to constraints 



and 



Z y. = 1.0 
1 



y . > 0.0. 

From previous discussion and assuming the p^ and 3^^ 
are known beforehand, this problem can be solved using the 
Direction Finding Algorithm. Observe that the first partial 
derivative required in the use of the algorithm is 



F (Yj = Z p. 3-- exp{-Z 3-*y-}. 

7i i iJ ^ j 13^3 

Solution of this problem via the Direction Finding Algo- 
rithm was accomplished on the IBM 360 computer system at the 
U. S. Naval Postgraduate School. Programming was done in 
FORTRAN IV Language using the G-Level Compiler. A copy 
of the program listing is included after the appendices. 

A. EXPLANATION OF EXAMPLE 

Since the purpose of this example is to demonstrate ap- 
plication of the algorithm and not to solve a given, existing 
problem, assignment of values to parameters of the problem 
was based on academic interest and convenience rather than 
practically oriented. 
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A grid size of 3600 squares, or equivalently a 60X60 
matrix was used. It is felt this grid sire selection is 
consistent with what might be used in a truly practical 
application and at the same time it enlarges the problem 
sufficiently so that information about machine time for 
calculation would be meaningful. 

The distribution of probabilities (p^'s) used in this 
example was uniform in nature. A p = 1/3600 was assigned 
to each of the 3600 grid squares. In a normal application 
it is expected that most grid squares would have a p = 0.0 
assigned to them and only a relatively small number would 
be assigned a positive probability. Those positive assign- 
ments would be most likely based on observer reports, knovm 
routes of movements, specific areas offering good cover and 
concealment, and other considerations such as these. How- 
ever, in this example a uniform distribution was selected 
for two reasons. Firstly, making all p’s positive is con- 
sidered a worst case condition with respect to the number 
of individual calculations required and the resultant ma- 
chine time necessary for execution. Since solution by the 
Direction Finding Algorithm is an iterative process, machine 
time is an important consideration and a known upper limit 
on this factor is extremely meaningful. Secondly, a uniform 
distribution was used so that effects on the boundaries of 
the grid could be observed. Since the damage function ' 

implies interaction between grid squares, it is not immedi- 
ately obvious what the optimal allocation along the edges of 
the grid should be. 
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Selection of an appropriate damage function would be 
based on the type of weapon being employed and the result 
of statistical data obtained for that weapon.. Inherently, 
the damage function selected directly affects total machine 
time necessary to solve the problem. However, this is a 
variable which is dependent on a specific situation and 
not of prime consideration in this paper. A simple expon- 
ential relationship was selected for the example because 
it was felt that an exponential function would be repre- 
sentative of most damage functions in terms of accuracy 
and machine time necessary for calculation.. 

Given the above explanation of problem parameters it 
now becomes necessary to develop a computer oriented logi- 
cal approach to solution of the problem. 

The logical sequence followed in finding the direction 
of fastest increase is explained earlier in this paper. At 
this point it should be noted that the program actually em- 
ploys two Direction Finding Algorithms. In one case 
(SUBROUTINE DF$ALL in the program listing) the y or direc- 
tion element associated with each and every element of the 
allocation vector is calculated. This is necessary before 
a final optimal solution can be obtained. However, while 
working toward the optimal solution, a significant reduction 
in calculation time can be realized by finding the direction 
of fastest increase only on the variables which have a posi- 
tive allocation currently associated with them or a positive 
Y in case the allocation is currently zero. For this reason 
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the program maintains a list (SUBROUTINE LIST) of those ele- 
ment indices which meet this criteria. With each iteration 
then, SUBROUTINE DF$LST is employed to find the direction 
of fastest increase only for those elements appearing on 
the "list.” When the optimal solution is found for the 
current "list," execution of the program shifts to SUB- 
ROUTINE DF$ALL. After execution of DF$ALL, a new^ list is 
formed and the above outlined process is repeated.. Only 
when optimization is indicated through the use of DF$ALL 
is execution terminated. 

In order to employ the Direction Finding Algorithm as 
indicated, various other calculations obviously must be 
made. Other routines in the program are written to accom- 
plish tasks such as evaluating the gradient, evaluating 
the objective function, and revising the current allocation 
to one which is closer to the optimal solution. These tasks 
must necessarily be carried out during each iteration of 
the solution process. A general flow chart outlining the 
logical procedure followed is depicted in Appendix A. A 
complete listing of subroutines employed along with a des- 
cription of the purpose of each is provided in Appendix B. 

Results obtained in the solution to this example prob- 
lem are included. It is interesting to note in these results 
that the optimal solution has no allocation of resources in 
the rows or columns immediately adjacent to the edges of the 
grid. Obviously the damage resulting in these grid squares 
from rounds allocated to and impacting in nearby interior 



22 



grid squares is considered sufficient to compensate for loss of 
damage effect that would occur if rounds were allocated to 
the outermost rows and/or columns. The Irigher allocation in 
those nearby interior grid squares suggests this conclusion. 

It is also interesting to note that the optimal solu- 
tion is not completely symmetric as might be expected from 
using a uniform type distribution of the p^*'s.. Inspection 
of the final allocation shows the upper left corner of the 
array somewhat different from allocations assigned in the 
other three corners. This is the result of the starting al- 
location used in the solution process and the precision con- 
sidered acceptable within the program. The initial starting 
feasible solution was for the entire inventory to be assigned 
to the extreme upper left grid square. This heavy allocation 
in just one location prevented a lesser positive allocation 
from appearing in nearby grid squares as the program stepped 
along towards optimality. The program was actually moving 
toward a totally symmetric final allocation, but settled for 
something less because of the precision cutoff programmed 
into the problem. If a higher degree of precision were de- 
manded and if the program were rewritten to provide greater 
machine accuracy, the program would tend to a completely 
symmetric optimal solution. 

Certain difficulties were encountered in programming 
this problem for machine solution. A discussion of some of 
these difficulties and comments on observations is considered 
of academic interest to the reader at this point. 



23 



i 




Since use o£ the Direction Finding Algorithm is an 
iterative process which moves from an.y feasible starting 
solution to a better one and then ultimately to an approxi- 
mate optimal solution, selection of an initial starting 
point can be a critical matter. Initial attempts at solving 
this problem were done by assigning the entire inventory 
to the first grid square and proceeding from there. Anal- 
ysis of results obtained from this "start from scratch" ap- 
proach indicated that machine time • in excess of twelve hours 
would be required to completely solve the problem. This re- 
quired time was considered unacceptable and, therefore, an 
alternate approach was devised. The sector covered with a 
3600 grid was initially considered to be covered with a 400 
element grid. An optimal solution was found for the 400 grid 
case and this solution was then used as a starting solution 
for the 3600 grid case. Results in terms of machine time 
were remarkable. The optimal solution to the 400 case was 
obtained in 17.5 seconds of machine execution time. Over- 
all machine time to the final optimal solution was six 
minutes nine seconds. 

A point of interest in the program is the value of "D" 
or the distance which the allocation is moved in each itera- 
tion. Recall the lower and upper boundaries of the alloca- 
tion elements are 0.0 and 1.0 respectively. The lower bound 
on D is 0.0 which is synonomous with no move at all. The 
upper bound on D is the maximum distance across the Y-space 
and is /~2. In some calculations within the program D is 
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controlled by a variable which is near its boundary and 
trying to move in the direction of its boundary.. The 
result of this situation in large problems such as this 
example is that D becomes very small and this leads to 
short moves on each iteration. When many "yrarlables are 
close to their boundaries and trying to move in that direc- 
tion, the overall impact is a significant increase in 
machine execution time. However, there are various places 
in the solution process \\rhere D can be successfully reset 
to a predetermined value so that time is not wasted on 
relatively small moves when larger moves are permitted. 

There is also the problem of what to do with D when move- 
ment is made in a direction of increase but it has not re- 
sulted in an improved value of the objective function. To 
resolve these difficulties, an initial D value of 0.5 was 
selected. This had the effect of making the first step a 
large one. A reset value of 0.1 was decided upon and each 
time a new list was compiled, D was reset to this value. 

A modification procedure was used to adjust D when a move 
was made but an improvement in the objective function was 
not realized (rightly enough referred to as a "failure”) . 

In this case D was simply cut in half and the move tried 
again. In the case of success in improving the value of the 
objective function, D was left alone and the same value used 
in the next iteration. The only exception to this procedure 
was in the case of a "success" which had been immediately 
preceded by a failure. In this circumstance D was halved 
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for the next iteration. To understand the rationale be- 



hind this consider the following simplified two-dimensional 
representation of a concave function which is to be maxi- 
mized : 




Let point A represent the current allocation and the direc- 
tion of fastest increase is obviously to the right. Sup- 
pose the present D is equal in value to (B-A) .. When the 
move to B is attempted, a "failure" occurs since the value 
of the function is less than it was at point A. The pro- 
cedure now calls for halving D and try to move again. This 
time the move is to point C and a "success" is realized. 

If D is not adjusted at this point, the next iteration 
would call for an attempted move back to point A which was 
the starting point for the previous iteration. To preclude 
this return to a previous allocation, D is cut in half 
after the occurrence of a success which had been immedi- 
ately preceded by a failure. It should be noted that this 
procedure becomes particularly useful when, the optimal 
solution is being neared and oscillations around the optimal 
solution start occuring. 
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Another consideration worthy of note is the influence 
of machine accuracy on the results of the problem. In a 
problem of this sort with potentially 3600 variables posi- 
tive simultaneously, the effects of calculating, founding, 
and cumulatively adding these variables must be taken into 
account. Early attempts at solution were made using single 
precision under FORTRAN IV. This approach was quickly 
changed to double precision for all real variables used in 

the program. Even so, results obtained were considered 

_ 1 0 

suspect in accuracy beyond 10 . To illustrate the prob- 

lem, in this example the initial p^ value assignment to 
each grid square was made by assigning each grid square a 

value of 1/400. This value is supposed to be accurate in 

_ m 

FORTRAN IV to 10 . As a check device the program then 

sums all p's assigned and prints the results. The re- 
sultant sum of 400 separate additions differs from 1.0 by 
”■ 8 

2.24x10 . Whether or not accuracy at this level is tol- 

erable is up to the user, but most certainly is worthy of 
some • consideration . 

Machine storage space and execution time requirements 
in this problem were sufficiently significant to be worthy 
of comment. Because of the requirement to maintain a large 
number of large arrays of stored information (in the common 
storage region alone there are six arrays each of which con 
tains 3600 double precision real numbers) , this example pro 
gram run on the IBM 360 System required 332,000 bytes of 
storage. Some space saving techniques were employed and. 
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undoubtedly, additional ones could be found. However, it 
should be recognized that the storage requirements are not 
a direct result of the methodology used in the Direction 
Finding Algorithm, but rather the result of the nature of 
the problem being solved. As mentioned previously, total 
execution time necessary was six minutes nine seconds. The 
approach of solving the problem on a small scale first and 
then using these results as a starting solution for the 
larger case was an instrumental factor in keeping execution 
time within reason. Other techniques involving the applica 
tion of "clean coding" contributed toward reducing time 
requirements. One technique developed is considered to be 
of academic interest and involved the calculation of geo- 
metrical distance between squares. This calculation is 
necessary in evaluation of the gradient and in evaluating 
the function value. The distance calculation involves sim- 
ply measuring the horizontal and vertical distance between 
the two grid squares in question and then evaluating the 
square root of the sum of squares. To illustrate the signi 
ficance of this simple routine, it is done in excess of one 
million times in each iteration of the program. The execu- 
tion time required to call the built-in machine square root 
function and evaluate the argument is known to be in excess 
of 200 microseconds. Since this obviously contributes to 
exorbitant machine execution time, it was decided to make 
the distance computation between any two grid squares once 
and for all at the beginning of the program and then store 
this information for later reference as required. This was 
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possible since the distance evaluations were somewhat re- 
petitive in nature. The exact time saving incurred by ap- 
plication of this technique was not evaluated. Hcnvever, a 
conservative estimate is that execution time was reduced 
by a factor of five and, therefore, makes the nominal ad- 
ditional storage requirement well worth while. The specific 
operation just described is contained in SUBROUTINE REFER. 
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IV. SUMMARY AND CONCLUSIONS 



The Direction Finding Algorithm can be applied to a 
variety o£ concave or convex differentiable fun_ctions to 
find an optimal solution. However, the domain of these 
functions must be represented as a closed convex, set. 

Its applicability in practical problems is particularly 
useful when the objective is to optimize allocation of 
resources. 

In large scale problems it is easily adaptable to 
machine utilization. Time and storage requirements for a 
machine solution are primarily dependent on the nature of 
the specific problem being solved. The algorithm has 
relatively little effect on machine storage requirements. 
However, the algorithm's impact on machine time is the re 
suit of the boundary problem previously discussed. A 
method of overcoming the boundary problem as well as an 
extended application of this algorithm into an. area of 
problems which require the allocated variable to be repre 
sented in typical matrix form rather than in vecfor form 
is currently being investigated by the aforementioned 
Dr. Danskin. 
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APPENDIX A: FLOWCHART -GENERAL LOGICAL PROCEDURE FOLLOWED 
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APPENDIX B: SUBROUTINE LISTING 



SUBROUTINE A>MCAL : AMMCAL evaluates the objective function. 

SUBROUTINE AREA: AREA determines which j 's are associated 

with a given i in the term. 

SUBROUTINE CHANGE; CHANGE converts the optimal solution 
obtained from the 400 grid case into the starting solution 
for the 3600 grid case. It also changes values of various 
parameters used within the program for use in the 3600 grid 
solution. 

SUBROUTINE CNTROL : CNTROL controls the sequencing of events 

in solution of the problem. It tests each attempted move 
for success or failure and adjusts distance (D) accordingly. 
SUBROUTINE DF$ALL: DF$ALL is the Direction Finding Algo- 

rithm which finds the direction of fastest increase for the 
entire allocation vector. 

SUBROUTINE DF$LST: DF$LST is the Direction Finding Algo- 

rithm which finds the direction of fastest increase just 
for those variables on the list. 

SUBROUTINE LIST: LIST compiles a list of those variables 

which are positive and/or those variables which have a 
positive y. 

SUBROUTINE MOVE: MOVE adjusts the allocation vector to 

reflect an attempted move in the direction of fastest in- 
crease. It moves to the boundary any variables which are 

.14 

within 10 of their boundary. 
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SUBROUTINE OMCALC : OMCALC calculates the vector of first 
partial derivatives of the allocation vector. 

SUBROUTINE REFER: REFER fills a matrix with values rep- 
resenting geometrical distance between two grid squares. 
SUBROUTINE SET$P: SET$P sets precision desired into the 
program and fills the P vector. 



33 



4 



COMPUTER RESULTS 
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THE INDICATES THE SUBSCRIPT NUMBER OF THE ALLOCATION VECTOR ELEMENT 

APPEARING IN THE INDICATED CORNER OF THAT SECTOR. 



BELOW IS LEFT SIDE OF OPTIMAL SOLUTION OF AOO GRID CASE 
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BFLnW IS RIGHT SIOE OF OPTIMAL SOLUTION OF AOO GRID CASE 
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NOW CONVERTING TO ^600 CASE 



BELOW IS SECTOR I OF OPTIMAL SOLUTION 
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Bf^LOW IS SBCTOR II OF OPTIMAL SOLUTION 
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BELOW IS SECTOR III OF OPTIMAL SOLUTION 
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cooccc ccccc COCCCCOCCCCCOCCOCCO 



ooc occ ccccc ccccc ccc cccoc coo ccc 

^ p-< P—4 P—4 1—4 P-4 p-4 ^ 1—4 p>4 p>4 fv CV CM vP sP ^P vP Vp >P 'P 'P O 

m (p. CP fp fp- fp fp fp. CP CP fp. fp, CP fp fp. CP. r' CP. nO vP 

occcccccccc ccccc ccccocco 
occ ccccc occccccc occccccc 
ccc CO ccc ccc ccccc occ occccccc occ 



OO O C coco ccc ccc C C ccc CO ccc ccoccc 

P—4 p— ’ ^ p-4 f— 4 p-4 p-4 p-4 p-4 p- 4 ^ f— (\| f\j s0 \Q \Q vO \C 'P O O 

fPfpifpiCPrPfPrp fpcp.fp.rprpcpfp rp fp rp (p^t<J■<^^P^P^P 

oooccccccocc cccc ccoccccc 
oco ccccc ccc ccccc ccc c-occc 
coc ccc ccc ccccc cc ccccc ooc cooccc 

oco occ ccc CO ccccc ccccc occ ccoccc 

^ ^-4 1—4 p-4 1—4 p-4 p-4 p-4 p-4 p-4 p-4 f\J (\J C\J \Q \Q \C vO 'P Vp ^ 

fPCPCPfp. CP rprp.rp fPfp.cPrPcp.cPrp>cPncPvJ-<t->J'>P'PvP 
ooc occccccc ccccc occccccc 
ooc ccc ccc CO ccccc occccccc 
occccccc CO cccc CO oc cccc cccoc ccc 



occccccc cccocc cc cccocc oc ccoccc 
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BELOW IS SECTOR XI OF OPTIMAL SOLUTION 



^>c^%£:^Ov0^r'^CsCvn>^^^•f^^*oooo'(TCr irir it' 

Lrir f^r-r- 

ocrcoocoocccccroccccccccoo 

ccoooooococoooooocccccco 

oococooocococooccococoococooco 



cooooooocccooooocooococcococoo 

%CNC'ONO>o%ooosc^-NO'£(^r^r^ooccr a cro^imr 

oooc ooococ c CC CCOCOOOCCOO' 
oooccc coocccocoococcoooo 
oocococcooccocccooocccocccoooc 

occoooooccoooocoocccocoocooooo 

<J>c^c•'<CvC>cl^lr^L^lrLrL^^r:'^c^oC7 cr cr o oc od coco 
fn r*‘ f<^ r*' r*' r ' r r' rr. r ' r<' LT IT Lp vc vC sO 

coocoocococcccococoooooc 
cccoccccccococcoccccccoo 
ocooococccccocococoocccoccoooo 



ooocooooococooooooooooococoooo 

sC^ccosCvoi^ir^LrLmrin>CNnvoacrC'ooo<rcccc 
rofT.(<' c^i c^, fT fT fT fPfo fT' LT IT ^ ^ ^ 

occccoccecocooccccococcc 

oooocococococoooccccococ 

coooccoccoccccoooocoooooooccoc 

ooc oooc cooooocooooo oooc ooooooo 

oOsCsCc^iT \r \ririr\r 0^0 Q o^Q^c: (DC CC CO CO 

CO r'. fp f^. r<^ TO rn r' r'. o ‘ f<"» r<^ fp Lp u" LT vO 

ocooocc ccccooccccoc eccoo 

coo oo c ocoooococ oooc coo oo 

coo ocoo oooc ooococ oo coco ooooooo 



oocooooocccocooooocoocoooooooo 

rvJCv^c\icv 4 CNj(\j(\icvojcv(\,(xjrororn nC vO >cr^r-r^xCNC vC 
f<^ fo fo, fo (T, rn r<^. sf sC’ vO sC 

coo coco oooc oooc c coo ocoo c 
ooooooo coco oooc CO ooooooo 
oococcocccccocooccooooooc oocoo 



coo oooc oo oo ooco ocoo coco occ ocoo 

cvjcvjc\jcv.’cvc\!cv(xjcvr0f\;cv:r<" rf' fT' >c sC 

rn o ro O'! r^. ro rn ro o, r'. fo ro o o r'. o vO vC 'O 
o ooooooo ocoo coo oocoo ocoo 
ooooooo ooc coo ocoo ocoo ooc 
OOOOcOC occ occ occ oocoo cc ooooooo 

ooo Co oooc ooooooo ooc ocoo ooooooo 

rs.TvcvrJcsjcvfvcNJcvrvjcvcv.r*' for'i vCvCvcr^f^r^sO'CiNC 

occ ooococ OOOOCOC ooc c c oo c 
ooo oooc ooc ooo coo ooooooo o 
ocooocc ocococooooccocoo ooooooo 

oooo ooococ ooooooo ooc ooo ooooooo 

C O O sO o 

ro ro r> rr (T'j (r^ fo rn CP rn ro rn o CO ro >0 nC >0 

ooococ cocccoocoococcccoc 
ooo oo oo o oc o oooo oooo ocoo c 
cccocccccccoooooccoococoooooco 

ooococ oooo ooo oooo ooc ooo ooooooo 

Jii--*C\lC\lC\|lPlPir'^0>CN0 vO >C nO 

cooocccooooooocc oooc oooc 
ococoocccooccccoooccooco 
oooo oo oo cc cc oooc ooc oooo oocoo oo 

oooo oooc ooco oooo oooo ooc ooooooo 
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BELOW IS SECTOR XII OF OPTIMAL SOLUTION 



cccccoceoccccccoooocccoocccocc 

oocccoccccoccoccooccocooccccoc 



oococccccooocccoooccccocooococ 

oocooococccocoocccoococococcco 



cooocccc cccccc ccocccococcccccc 
ooooocc ococcoccccccoococccooco 



ooococccooccoccooococococccoco 

oooccoocccooocoococccooaccccoc 



ccococo ccccocccocooocooococcco 
C OOC OCO CO CC O C CO coco cc cccccc coo 



cocccccccocc cccoccc cccco cccccc 



coccccc ccccccccoc cccccc ccccccc 

^>OsOvCNONOsCNC>CsC'C>CvC>c^ccoanccirLrirr^r''r^ 
>C'C'£)vC>c>o>i:vC'C'CNCCsc>o>Csc>£:vcr^i^r^oocx:cc 
ccccccccccccccccccccccco 
ccccccc ccccocccccc cccccc 
cccccccccccccccccoccccc ccccccc 

ccccccoccccccccoccccocccccccoc 

^r'CsO'*CsCvONOvCvO>CNC^CsCNONOooa:jocLririrr*^l^r^ 

cocccccccccocccccc ccc ccc 
cccccc cccccc cc cccccc ccc c 
cocccccccocc cccccc c ccccocccccc 



ccccccc cc ccc ccccocccccc ccccccc 

^^CvC'£>OsCsCso>c>o>ONCNCNOocooocirLpirr^r-f^ 

^s0n0^nCsCvC>c^ c vCvOnOnCc nCnC vor- r-r^ccocoo 
CCCCCC cccccc ccc ccccccccc 
cccccc c cc ccccccc cccccc cc 
cc cccccc ccccccccccoccc ccc ccccc 

cccccc ccc ccccccc ccccccccc ccc cc 

<?o<)%o<?>C'CvC'0'C>o>ch-f^r-occc‘C'CrLrLnu:i' 

ccccccc ccccc ccc ccccccccc 
ccccccc ccccc cocccc cccocc 
cccccc cccccc cccccccc ccc ccccccc 

coco cccccc cccccc ccc ccccc cccccc 
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HEPE IS THE MAIN PR3GRAM 
YY IS THE BASE dOINT 

J IS A POINTER TO ELEMENTS ON THE LIST 

LONG IS THE TOTAL NJMBER OF PLEMENTS IN VECTORS 

SIZE IS LENGTH OF ONE SIDE OE SQUARE GRID 

EAILI IS TRUE IF SUCCESS COLLOWS A FAILURE 

SHORTD IS TRUE IP WORKINO WITH A D CLOSE TO THE BOUNDARY 

jjj IS 1 IE abbreviated omcalc is desired 

IMPLICIT REALM'S! A-H,o-Z,$) 

LOGICAL SHORTD, FATLl 

■ COMMON YI3600) , VV ( B 600 ), J ( 3600 ) ,GA(3600) rOM(3600) , 
1IM(3600) ,P(3600) ,E( 3600 K DSTNC E ( 8 , 8 ) , D , V? , AMM , AM, KO , 
ILONGtll I ,JJJ, JKtISIZP, I RANGE, ILO, I MI , i T ER , SHORTD, E A IL 1 
PLANK COMMON STATEMENTS IN SUBROUTINES HAVE BEEN OMITTED 
FOR THIS printing 
11 = 0 
T 11 = 0 
JJJ=0 
KKK=0 
KKKK=0 
NN=0 
ITER=0 
L0NG=400 
ISIZE=20 
IPANGE=2 
IL0=3 
IHI=18 
AM=1.0 
D=0.1 

FAIL1=. FALSE. 

SHOPTD=. FALSE. 

DO 1 1=1,3600 
YU ) =0.0 
YV( n=n.o 

OMU )=0.0 
GA( I )=0.0 
TM( I )=0 
J( I )=0 

1 CONTINUE 

FOLLOWING TWO CAROS INDICATE STARTING SOLUTION 
Y( 1)=1.0 
YYm = 1.0 
CALL SPT^P 
CALL REFER 
CALL CNTROL 
CALL C.MANGE 
CALL CA'TPOL 

AT THIS POINT problem IS SOLVED Z NECESSARY WRITE 
STATEMENTS NEED TO BE INSERTED 
STOP 
END 



SUBROUTINE CNTROL 

THIS PORTION CONTROLS EXECUTION OF THE PROGRAM 
IMPLICIT RPAL^=8( A-H,0-Z, $) 

LOGICAL SHOPTOtFAILl 
1000 CALL LIST 

CALL AVMCAL 

FIRST TIMP THROUGH AMMCAL IT IS CUT OPE SMQRT 
CALL Of-' CALC 

CALL OP5ALL(Y,OM,GA,LONG,V2, S1090) 

AT THIS POINT GET ON WITH NORMAL EXECUTION 
1010 CALL LIST 
0 = 0.1 

IF(III.EO.O) 3=0.8 
I 11=1 

FAIL1=. FALSE. 

SHORTD=. FALSE. 

1015 DO 1020 K=1,K0 
JK=J(K) 

YYY=YY( JK) 
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GAA=GA( JK) 

IF(YYY+n*GAA.GE.O.O) GO TO 1020 

D=-YYY/GAA 

SHOPTD=.TRUE. 

1020 CONTINUE 
1025 CALL ‘■invE 

CALL ArMC^L 

IF(A»"M.LT,AM) GO TO 1040 

BRANCH TO 104v0 BELOW IP SUCCESS-OTHEPVfl SE DO FOLLOWING: 

. 0 = 0/2 
FAIL1=. TRUF. 

GO TO 1025 

IF SUCCESStPI'^ST SET NEW AM C HOVE BASE POINT 
1040 AM=AMM 

DO 1045 K=1,K0 
JK=J(KI 

1045 YY(JK)=Y(JK) 

IF(FAILl) GO TO 1050 
GO TO 1055 
1050 n=D/2 

FAIL1=. FALSE. 

1055 CALL OMCALC 

CALL DF$LST( V, OM,GA,Jt KO,LONG, V2,C1060) 

1F{ .NOT.SHORTD) GO TO 1015 

0 = 0.1 

SHORTD=.FALSF. 

GO TO 1015 
1060 JJJ=1 

here calculate only those omegas not ALREADY CALCULATED 
DO 1070 JJ=ltLONG 
IF( IM( JJ I.EQ.l ) GO TO 1070 
0M( JJ)=0.0 
JK=JJ 

CALL OMCALC 
1070 CONTINUE 
JJJ=0 

CALL 0F$ALL( Y»OM,GA,LONG,V2,S1090) 

GO TO 1010 
1090 RFTURN 
END 



SUBROUTINE SFT5D 

THIS ROUTINE ENTERS P VALUES AT START OF PROGRAM 
IMPLICIT R'^AL*8( A-H,n-7, $) 

LOGICAL SHORTDtFAILl 
PREC=10.0=^-(-3 ) 

V=PPEC/2=<'*. 5 
V2=V=!=V 

20 • VAL=1. 0/400.0 

DO 22 N=l,400 
22 P(N)=VAL 
SUMP=0. 0 
DO 25 N=l,400 
25 SUMP=SUMP+P(N) 

DO 27 N=l,400 
27 P(N)=P( NI/SU^P 

RFTURN 
END 



SUBROUTINE LIST 

THIS IS LIST CONSTRUCTING SUBROUTINE 
IMPLICIT RFAL*B( A-H,n-Z,S) 
LOGICAL SH0RTD,FAIL1 
K=0 

ITER=ITEP+1 
DO 60 JJ=ltLCNG 
IF(Y( JJ I .GT.O. 0) GO TO 50 
IF(GA( JJ) .GT.O ,0) GO TO 50 
GO TO 50 
50 K=K+1 
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J(K)=JJ 
60 CONTINUF 
KO=K 

KO NOW CONTAINS LIST SIZE 
0 = 0.1 
RETURN 
END 



SUBROUTINE DF$L ST {$, $, r,A$ , J$ t K0$ , LONGS » V2$', * ) 

THIS IS DIRECTION FINDING ALGORITHM USING THE LIST 
OFSLST ACCEPTS AS paRANETEoS THE FOLLOWING: 

S IS VECTOR OF VARIABLE TO BE ALLOCATED 

OMS IS VECTOR OP FIRST PARTTALS 

GAS IS GAMMA VECTOP 

LONGS IS CURRENT VECTOR LENGTH 

V2S IS V*V WHERE V= PR FC I SI ON / 2 . E 

JS IS VECTOR OF POINTERS tq ELEMENTS ON THE LIST 

KOS IS A number indicating the LIST LENGTH 

=<' IS SPECIAL RETURN POINT (RETURN 1) 

IMPLICIT RPAL^=R( A-HtO-Z,S) 

REAL*8 LA, MU 
LOGICAL IND 

DIMENSION S(LONGS) , CM S ( LONGS V, GA $( LONG $ ) , JS(LONOS) 
THIS SUBPOUTINE P<=9UIRES NO BLANK COMMON ACCESS 
C0MM0N/W0RKA/IMD(B600) 

100 N=0 

KK=0 
S = 0.0 

DO 110 K=1,K0S 
JK=JS( K) 

GA$( JK) =0.0 
IND(JK) =. FALSE. 

IF( S{ JK ) .EO.0.0 ) GO TO 110 
IF(5(j:<).E0.1 .0) GO tq 110 
INO(JK)=.TRUE. 

N=N+1 

S=S+nMS ( JK) 

110 CONTINUE 

IF(N.EQ.O) GO TO 115 

MU=S/N 

AOM^rMU 

GO TO I AO 

115 AOM=PM$ ( JS( 1 ) ) 

DO 120 K=1,K0$ 

JK=JS( K) 

OMM=OMS( JK) 

IFCOMM.LT. AOM) AOM=OMM 
120 CONTINUE 
MU=AOM 
lAO JJ=0 

DO 150 K=1,K0S 
JK=JS(K) 

IFdND(JK)) GO TO 1^0 
IF( S( JK) .LT.l. 0) GO TO 150 
1F(0MS( JK) .GE. AOM) GO TO 150 
AOM=OMS ( JK) 

JJ = JK 

150 CONTINUF 

IE( AOM.GE.MU) GO TO 170 
KK=0 

MU=(N*MU+A0M)/(N+1) 

AOM=MU 

N=N+1 

IND( JJ)=.TRUE. 

GO TO lAO 
160 JJ=0 

DO 165 K=1,KCS 
JK=JS(K) 

IE( IND( JK) ) GO TO 165 

IE( $( JK) .GT.0.0) GO TO 165 

IF(DM$( JK) .LF. AOM) GO tq 165 
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AOM=OM${ JK) 

JJ=JK 

165 CONTINUE 

1F( AOM.LE.MU) GO TO 130 
KK=0 

MU=(N«"^U+AOM)/ ( N+1 ) 

AOM=MU 
N=N + 1 

INO( JJ)=.truE. 

GO TO 160 

170 IF(KK.EQ.l) GO to 190 
GO TO 160 
180 KK=1 

GO TO 140 
IPO SS=0.0 

00 195 K=1»K0$ 

JK=J$(K) 

IF( .NOT. INDUO ) GO TO 195 
Sn=OM$( JK)-NU 

ss=ss+sn*sD 

195 CONTINUE 

IP(SS.LT.V2$) RETURN 1 
LA=DSORT(SS ) 

DO 200 K=1,K0$ 

JK=J$(K) 

IF( .NOT. INDUO ) GO TO 200 
SD=0M$( JK)-MU 
GASUKI =SD/LA 
200 CONTINUE 
RETURN 
END 



SUBROUTINE PFSALU « t99 5,r,A$,L0NG5,V?4,T) 

DF$ALL ACCEPTS AS PARA^FTEPS THE FOLLOWING: 

$ IS VECTOR OF THc VARIABLE tq BE ALLOCATED 

OMi IS VECTOR OF FIPST ?AC»TTALS 

GA-5 IS GA^^MA vector 

LONGS IS CURRENT VECTOR LENGTH 

V2S IS V*V WHERE V=PREC I STON/2*^. 5 

* IS SOECTAL Return POINT ( REtupn 1) 

IMPLICIT PEALTSI A-H,0-Z,$) 

RFAL’i'S LA, MU 
LOGICAL TND 

dimension S (LONGS) , OM S ( LONGS ), G A S ( LONGS ) 

THIS SUBROUTINE REQUIRES NO BLANK COMMON ACCESS 
COMM0N/WORKA/INn(3600) 

IND(I) IS SET tbuF WHEN S(I) ENTERS DS 
500 N=0 

N IS COUNTFR FOR NUMBER OF ELEMENTS IN DS 
KK=0 

KK IS LOGIC SWITCH WHICH CONTROLS FLOW 
S=0.0 

FIRST TO CHECK FOP ELEMENTS MOT ON THEIR BOUNDARIES 
DO 510 JJ=1, LONGS 
GAS( J J) =0. 0 
IND(JJ)=. FALSE. 

IF( $UJ ) .EO.O. D) GO TO 510 
IF(S( JJ ) .EO.l.O) GO TO 510 
IND( J J)=.TRUE. 

N=N+1 

S = S + OMSUJ) 

510 CONTINUE 

IF(N.EO.O) GO TO 515 
MU=S/N 

AT THIS STAGF, MU IS AVERAGE OF INTERIOR ELEMENTS 
AOH=MU 
GO TO 540 

IF DS IS EMPTY set-find MINIMUM OMEGA 
515 AOM=OMS( 1 ) 

DO 520 JJ=1, LONGS 
OMM=OMS{ JJ) 
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IP(OMM.LT.AOM) aom=omm 
520 CONTINUE 
MU=AOM 

NOW TO CHECK FOR ELEMENTS ON U°PER BOUNDARY AND NOT IN DS 
5A0 JMARK=0 

DO 550 JJ=1,L0NC$ 

IF( IND( JJ) ) GD TO 550 
IF( $( JJI .LT.1.0) GO TO 550 
IF(0M5{ JJ) .GE. AOMI GO tq 550 
AOM=OM$(JJ) 

JMARK=JJ 
550 CONTINUE 

IF( AOM.GE.MU) GO TO 570 
KK=0 

NOW TO pfcALCULATF MU 

MU=(N=i'MU+AOM)/ (N+1 ) 

AOM--MU 

N=N+1 

IND( JMARK)=.TRUE. 

GO TO 540 

NOW TO CHECK fqP ELEMENTS ON LOWER BOUNDARY AND NOT IN DS 
560 JMARK=0 

DO 565 JJ=l, LONGS 
IF( IND( JJ) ) G3 TO 565 
IF{S<JJ).GT.O.O) GO TO 565 
IF{OMS( JJ) .LE.AOM) GO TO 565 
AOM=OMS(JJ) 

JMA»K=JJ 
56 5 CONTI N'.IF 

IF( AOM.LE.MU) GO TO 5R0 
KK=0 

RECALCULATE MU AGAIN 

MU=(N*mu+AOM)/ (N+1) 

A0M=MU 

N=N+1 

IND( JMARK)=.truE. 

GO TO 560 

570 It=(KK.F0.1) GO TO 5Q0 

KK=1 IMPLIES NO MORE SEARCHING IS PEOUTRED 
GO TO 560 
580 KK=1 

GO TO 540 

NOW TO PINO SUM OP SQUARES OP ELEMENTS IN DS 
590 SS=0.0 

DO 595 JJ=1, LONGS 
IF( .NOT.INOI JJ ) ) GO TO 595 
SD=OM$( JJ)-MU 
SS=SS+SD*SD 
595 CONTINUE 

IF(SS.LE.V2S) RETURN 1 
LA=DSOPT(SS ) 

LAMBDA EQUALS SQUARE ROOT OF SUM OF SQUARES 
DO 6C0 JJ=1, LONGS 
IF( .MOT. IND( JJ ) ) GO tq 600 

$(I) NOT IN DS IMPLIES GAMMA(I)=0 
SD=OMS( JJ )-MIJ 
GAS( JJ) =SO/LA 

GAMMA VFC.TOR IS FILLED HERE 
600 CONTINUE 
•RETURN 
END 



SUBROUTINE AMMcaL 

THIS PART CALCULATPS THE H'S AND AMM 
IMPLICIT RPAL*8 ( A-H,0-Z, $) 
LOGICAL SHORTO,FAIL1 
DIMENSION H(3600) ,0(3600) 
COMMON/WOR'<B/H ( 3600 ) , 0 ( 3600) 
COMMON/ XFER 1 /L T M,LL ,LM, ML, MM 
300 DO 305 1=1, LONG 
305 H(I)=0.0 
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DO 370 K=1tK0 
JK=J(K) 

IF{Y{ JK) .EO.0.0) GO to 370 

CALL AREA 

DO 360 LI=LLtLY 

DO 350 MT=MLtMM 

I = (L1-1 )<=IST7E+MI 

IDL=TABS(LT-L) 

IDM=1 ARS(MI-M) 

IDLS=inL+l 
1DM$=IDM+1 
S=DSTNCE( 10L$» TDM$I 
B=DEXP(-S) 

1F{LONG.EQ.3600) GO TO 341 
1F(S.GT.3.5) B=0.0 
GO TO 342 

341 IF(S.GT.10.5) B=0.0 

342 H(! >=Hf T)+B*Y( JK) 

350 CONTINUE 

360 CONTINUE 
370 CONTIMUF 

DO 375 1=1, LONG 

1F( P{ I ) .‘"Q.0.0 ) GO TO 374 

Q{I)=DEXD(-Hn )) 

E(I) = Pm=J=0(I) 

GO TO 375 
374 E{I)=0.0 

37*^ CONTINUE 

IF(III.EO.O) RETURN 
S=0.0 

DO 3B0 1=1, LONG 
380 S=S+E(1) 

AMM=S 

RETURN 

END 



SUBROUTINE OMCA?.C 

THIS ROUTINE OALCULATES thE ONEGAS 
IMPLICIT R<=AL4«( A-H,n-Z,$) 
LOGICAL SH0PT0,EAIL1 
COMMON/ XFER 1 /L , M, LL , LM , ML , MM 

400 IFIJJJ.EO.l) GO TO 420 
DO 410 1=1, LONG 
0M(I)=0.0 

410 IH{I)=0 

IF( IIl.NF.O) GO TO 415 

DO 490 JJ=1,LGNG 

K=K0 

JK=JJ 

GO TO 420 

415 DO 490 K=1,K0 
JJ=LCNG 
JK=J(K) 

420 CALL AREA 
IM{ JK)=1 
DO 480 LI=LL,LM 
DO 470 MI=ML,mm 
I=( L l-1 )^ISIZE+MI 
IF(P( I ) .EO.O.O)GO TO 470 
inL=I ARS(LI-L) 

IDM=IAPS(MI-M) 

IDL$=1DL+1 
1DM$=IDM+1 
S=DSTNCEUDLS, IOM$) 

B=DEXD( -S) 

l«=(LrNG.E0.3600) GO TO 441 
IFIS.GT.^.5) ?=0.0 
GO TO 442 

441 IPf S.GT.10.5) B=0.0 

442 0M( JK)=CM(JK)+B«E{T) 

470 CONTINUE 
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480 CONTINUE 

IF( JJJ.EQ.DRETURN 
490 CONTINUE 
RETURN 
END 



810 

815 



820 

825 



830 

835 



840 

845 



SUBROUTINE area 
IMPLICIT REAL«8( A-HtO-Z,$) 
LOGICAL SHODTDfPAILl 
COMMON/ XFFR 1 /L »M,LL tLMt ML, MM 
L=( JK-1 ) /ISIZE+1 
M=JK-IS1ZE*(L-1 ) 

TF(L.LT.ILO) 30 TO BIO 

LL=L-IPANGE 

GO TO 315 

LL=1 

IFIL.GT.IHT) GO TO 820 
LM=L+I RANGE 
GO TO 825 
LM=ISI?F 

IF(M.LT.ILO) go to 830 
ML=M-I RANGE 
GO TO 835 . 

ML=1 

1F(M,GT.IHT) GO TO 840 

MM=M+IPANGE 

GO TO 845 

MM=ISIZE 

RETURN 

END 



SUBROUTINE ^cv? 

IMPLICIT RPAL*8( A-H,n-Z,« ) 

LOGICAL SHORTOtFAILl 

30 DO 35 K=1,K0 
JK=J(K) 

YYY=YY( JK) 

GAA=GAUK) 

Y{ JK)=YYY+0*GAA 

IF(Y( JK.) .GT. 1.00-14) GO TQ 34 

Y{ JK)=0.0 

GO TO 35 

34 IF(Y( JK).LT.,PQ99P99999Q999) GO TO 35 
Y( JK)=1.0 

35 CONTINUE 
RETURN 
END 



SUBROUTINE REFER 

THIS ROUTINE LOADS OSTNCP MATRIX 
IMPLICIT RFAL*B( A-H,0-Z,S) 
LOGICAL SH0RT0,FAIL1 
70 DO 75 1=1,8 
DO 74 M=l,8 
T$=I-1 
M$=M-1 

-SS = I $+M$'l=Mt 

S=nSORT{ SS) 

74 OSTNCE( T,M)=S 

75 CONTINUE 
RETURN 
END 
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SUBROUTINE CHANGE 

THIS ROUTINE CONVERTS noTIM&L SOLUTION FROM ^00 CASE TO 

STARTING SOLUTION FDR 3600 CASE 
IMPLICIT REAL^^-BC A-H,0-Z,$) 

INTEGFR*4 BIGPOW,BIGCOL 
LOGICAL SHORTD^FAILl 
C0MM0N/W0RKB/PP(3600) , FILL (3600) 

D=0. 1 

PREC=10.0=«'*(-5) 

V=PPEC/2=^==!=.5 

V2=V=i'V 

AM=1.0 

L0NG=3600 

1ST7.E = 60 

IRANGE=7 

ILO=B 

IHI=53 

ITER=0 

II 1=0 

SH0RTD=. FALSE. 

FAIL1=. FALSE, 

DO 13 1=1,3600 
P( I ) = 1. 0/3600.0 
0M( I )=0.0 
GA(I)=0.0 
IM( I )=0 
J( I )=0 

13 CONTINUE 

DO 5010 1=1,400 

FIRST DETFOMINE ROW AND COLUMN OF INDEX NUMBER 
TFMPY=Y(! ) /o.o 
NRnw=(I-l)/?0+l 
NCOL=1-20* ( NPQW-l ) 

NOW TO FIND BASE SQUAPF WHICH IS TOP LEFT SQUARE 
BIGR0W=NRnw=?=3-3 
RIGCGL=NC0L43-? 

IBASF=PIGC0L+60*(BIGR0W-l) 

YY(IBASE) =TpMPv 
YY(IBASF+1) =TPMPY 
YY(TPASE+2) =TFmpy 
YY(IBASF+ 60) =TFMPY 
YYnBASE+61) =TFMDY 
YY(IBASP+62) =TFYdy 
YY( TBASE+120) = TE*^PY 
YY( IBASE+121 )=TFMPY 
YY( IBASF+122 )=TFMPY 
5010 CONTIN'UE 

NOW TO fill Y S P MATRIX 
DO 5020 1=1,3600 
5020 Y( I )=YY( I ) 

RETURN 

END 
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